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Dwarf galaxies pose significant challenges for cosmological models. In par- 
ticular, current models predict a dark matter density that is divergent at the 
center, in sharp contrast with observations which indicate an approximately 
constant central density core. Energy feedback, from supernova explosions 
and stellar winds, has been proposed as a major factor shaping the evolution 
of dwarf galaxies. We present detailed cosmological simulations with sufficient 
resolution both to model the relevant physical processes and to directly assess 
the impact of stellar feedback on observable properties of dwarf galaxies. We 
show that feedback drives large-scale, bulk motion of the interstellar gas re- 
sulting in significant gravitational potential fluctuations and a consequent re- 
duction in the central matter density, bringing the theoretical predictions in 
agreement with observations. 

Dwarf galaxies are the most common galaxy type 0. In the hierarchical picture of cos- 
mic structure formation, dwarf galaxies form first, later becoming building blocks for larger 
galaxies. Thanks to their proximity in the local universe (around 18 galaxies are located within 
300 kpc of the Sun), several of these galaxies have been studied in detail. By measuring line-of- 
sight velocities for hundreds of stars, accurate modeling of the mass distribution has revealed 
features that pose severe challenges for the standard cosmological model. It appears, for ex- 
ample, that the distribution of dark matter (which is the dominant mass component of these 
galaxies) is of almost constant density in a central region that is comparable in size to the stellar 
body of the galaxy (|2]|2]II]|5]|3|). In the best studied systems, Fornax and Ursa Minor, the radius 
of this region is ~ 400 pc and ~ 300 pc, respectively (6). This core is at odds with existing cos- 
mological models, which reliably predict the dark matter to have a divergent density (a cusp) at 
the galactic center ©. Some dwarf spheroidal galaxies also exhibit radial gradients in the stellar 
population, with stars more deficient in heavy elements (and therefore presumed older) having 
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a more extended distribution and being kinematically warmer than more metal-rich stars d#]|9]). 
Further, the presence of globular clusters in many dwarfs is puzzling since these massive, com- 
pact systems of many thousands of stars would have suffered gravitational drag as they moved 
through the dark matter background of the galaxy halo. This dynamical friction would have 
caused the globular clusters to spiral into the galaxy center on time scales much shorter than the 
age of the galaxy (the Fornax dwarf spheroidal galaxy is notable in this regard ©). 

It is well established that massive stars inject large amounts of energy into the surround- 
ing medium via stellar winds and supernova explosions, resulting in large-scale (hundreds of 
parsecs in large galaxies) random bulk motions of the interstellar gas at close to sonic speeds 
(~ 10 km s~ 4 for the typical gas temperature of 10 4 K) (|70][77II72|) . The effect of such per- 
turbations is larger for dwarf galaxies, as they have lower gas pressure due to the lesser depth 
of their gravitational potential wells. This stellar feedback has been invoked to explain at least 
some of the puzzling properties of dwarf galaxies. In particular, there has been considerable 
debate as to whether or not it can turn the theoretically predicted central dark matter cusp into a 
core G3\M\n\\M. 

Previous theoretical work has included both non-cosmological and cosmological model- 
ing. High resolution, non-cosmological numerical models with detailed descriptions of relevant 
physical processes (TiOl 1771 17^1) suffer from unrealistically symmetric initial conditions and a 
static description of the dark matter potential, and from the lack of gas accretion from the ambi- 
ent cosmic medium. Prior attempts at self-consistent hydrodynamic cosmological simulations 
have tended to focus on the formation of the very first small galaxy progenitors (TJ9ll20|) . or 
on dwarf galaxy models without sufficient resolution or the relevant physics to properly model 
star formation and feedback because of the substantial computational challenges involved in 
self-consistent modelling (|2ill22l) . 

Here we present the results of cosmological simulations of dwarf galaxy formation and evo- 
lution that adequately resolve and model the processes of star formation and stellar feedback. 
In good agreement with our previous semi-analytical results (1761) . our self-consistent model 
demonstrates that in small galaxies random bulk gas motions driven by stellar feedback play a 
critical role in determining the structure of the galactic center. The key result is the transforma- 
tion of the central density profile from a cusp to a large core. This is a consequence of resonant 
heating of dark matter in the fluctuating potential that results from the bulk gas motions. We 
also demonstrate that the same mechanism can explain other puzzling features of dwarf galax- 
ies, such as the stellar population gradients, low decay rate for globular cluster orbits and the 
low central stellar density. 

The simulations were run using the cosmological parallel tree code Gasoline (1231) . This 
code represents dark and stellar matter as a collection of dark matter and star particles and uses 
the smoothed particle hydrodynamics (SPH) formalism to describe gas evolution. A detailed 
description of the code, including the prescriptions for star formation and supernova feedback, 
can be found in refs. (1241 1251) . The very high resolution achieved in our models required the ad- 
dition of two key features to the standard cosmological code. First, low temperature (< 10 4 K) 
radiative cooling from de-excitation of fine structure and metastable lines of heavy elements 
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Figure 1: A zoom-in on the central part of the simulated forming dwarf galaxy at redshift 
z = 5.3. This time was chosen to illustrate the very clumpy gas distribution following a star 
burst. Gas is shown in blue and stars in red. (A) Global view. (B) View of the galaxy. (C) The 
central part of the galaxy. A number of star clusters are visible in panel C, the oldest (marked) 
has an age of 200 Myr. 

was necessary to correctly model gas cooling in small galaxies (f25l) . Second, since our mass 
resolution (< 200 M ) is sufficient to resolve individual supernovae, we introduced a new, 
stochastic, prescription for stellar feedback (|25l) . 

We created cosmological initial conditions with input constraints designed to produce a 
dwarf galaxy with total mass ~ 10 9 M at redshift z = 6 within a box of size 4 co-moving 
Mpc. A central, high resolution, sphere with radius 0.4 co-moving Mpc was populated with gas 
particles. The particle masses inside the high resolution sphere were 1,900 M Q for dark matter 
and 370 M Q for gas. The mass of particles generated to represent stars was ~ 120 M . At the 
end of the simulations, the total numbers of dark matter, gas, and star particles were 1.1 x 10 7 , 
4.5 x 10 6 , and 4.5 x 10 5 , respectively. The gravitational softening length was held constant at 
12 pc. Further model details, including the description of numerical convergence tests and free 
parameter studies, can be found in (|25l) . 

Two primary simulations were run. The first one included all the key physical effects: gas 
dynamics, star formation, and stellar feedback. (This was by far the most computationally 
expensive, consuming 6 x 10 5 CPU hours.) The second one was a dark-matter-only control 
simulation. The simulations started at z = 150, and ended at z = 5. 

In the simulations, the matter distribution develops the classic web-like or filamentary struc- 
ture on large scales, with the most massive galaxy forming around z = 10 at the intersection 
of the major filaments near the center of the computational box (Figure 1). The evolution of 
the galaxy is relatively smooth (no major mergers) between z = 8 and 5. The star forma- 
tion in the dwarf galaxy is very bursty, with major star bursts repeating approximately every 
~ 80 Myr, consistent with the non-cosmological models of (|70|) . Stars form predominantly 
in clusters, but many of them quickly disperse. Starting at z = 6.2, when the galactic stellar 
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mass reaches ~ 10 7 M , clusters which survive until the end of the simulation start forming. 
These long-lived clusters have broadly the same sizes (~ 10 pc, essentially unresolved in our 
simulations), masses (~ 10 5 M ), and heavy element abundance (~ 3% of solar) as globular 
clusters observed in the local universe. It is noteworthy that in the Local Group no old (early- 
type) dwarfs with stellar mass < 10 7 M have globular clusters, whereas all brighter dwarfs 
(with the exception of M32, which is severely tidally stripped by its host galaxy, M31) have 
globular clusters ©. This suggests that a galaxy has to be large enough (> 10 7 M in baryons) 
to produce globular clusters, in good agreement with our simulations. 

Feedback from the bursty and clustered star formation results in a dramatically perturbed 
interstellar gas distribution on large scales (hundreds of parsecs; Figures 1 and 2 and Movie SI). 
This is consistent with the observed (irregular) distribution of gas in dwarf galaxies ©. Note 
that at high-redshift this feedback does not expel gas from the galaxy, in contrast to the maxi- 
mum stellar feedback mechanism (|74l) . Instead, supernova explosions compress gas into large 
shells and filaments, which are confined to the central part of the galaxy and move with speeds 
~ 10 — 20 km s _1 (comparable to the speeds of dark matter particles). We showed previ- 
ously (|731) that gas motion with these characteristics results in efficient gravitational heating of 
the central dark matter and flattening of the cusp. 

The heating of dark matter in our model dwarf galaxy is highly effective (Figures 2-4). 
Whereas both density, p, and velocity dispersion, a, of the particles are strongly affected by 
the variable content of gas and stars at the galactic center, the phase- space density, F = p/cr 3 , 
is much less sensitive to adiabatic compression of dark matter by baryons (Figure 2). In the 
dark-matter-only simulation, F stays approximately constant, whereas in the hydrodynamic 
simulation, F gradually decreases with time as the result of the stellar feedback, becoming a 
factor of 10 lower than for the dark-matter-only simulations at the end of the evolution. 

The dark matter density is strongly affected by the stellar feedback only in the central region 
of the galaxy (Figure 3). This is the region in which the enclosed gas mass occasionally domi- 
nates that of the dark matter, and is where the gas is most strongly affected by the feedback. At 
the end of the hydrodynamic simulations, the dark matter density at the smallest resolved radius 
becomes a factor of seven smaller than in the dark-matter-only simulations. 

Whereas the dwarf galaxy halo in the dark-matter-only simulation develops a central cusp 
with logarithmic slope of 7 = —0.95, consistent with previous predictions of the standard 
model ©, in the hydrodynamic simulations, resonant heating due to stellar feedback turns the 
cusp into a flat core with radius 400 pc (Figure 4) and average density 0.2 M pc~ 3 . These 
core parameters are close to those inferred for Fornax, ~ 400 pc (6) and ~ 0.1 M pc -3 ©, 
respectively. The same mechanism produces a core of somewhat smaller radius (~ 300 pc) in 
the distribution of stars, and, significantly, pushes newly formed globular clusters away from the 
galactic center. The four oldest globular clusters, for example, were born with radial distance 
dispersion a r = 37 pc (essentially at the galactic center), but after ~ 200 Myr of evolution this 
distance had grown to a time-averaged value of a r = 280 pc (comparable to the stellar core ra- 
dius). We suggest that resonant gravitational heating can at least partially explain why globular 
clusters in Fornax, and in some other dwarfs, are located at large distances from the galactic 
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Figure 2: Evolution of the central quantities in the model dwarf galaxy. In the upper panel, 
solid lines correspond to changes in the dark matter (black), gas (blue), and stellar (red) masses 
enclosed within the central 100 pc as a function of the redshift, z. The dashed blue line shows 
the evolution of the enclosed gas mass within the central 1.6 kpc (half the virial radius). In 
the lower panel, green and black lines show the evolution of the central dark matter phase- 
space density, F, for the hydrodynamic and dark-matter-only simulations, respectively. We 
also show the evolution of the velocity anisotropy, 77, for the same dark matter particles as 
were used to calculate F (magenta line; horizontal black line marks 77 = 0). Here 77 = (of — 
a t)/( a r + °f)' wner e u r and a t are, respectively, the one-dimensional radial and tangential 
velocity dispersions. 
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Figure 3: Evolution of the enclosed dark matter masses in the model galaxy at different radii. 
Dashed lines correspond to the dark-matter-only simulation, and solid lines correspond to the 
hydrodynamic simulation. 
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Figure 4: Radial profiles for the model galaxy at redshift z = 5.2. At this time the central 
gas density is very low, minimizing the adiabatic compression of dark matter due to baryons 
(which makes it appropriate for comparison with presently observed gas-poor dwarfs). Green 
and red lines show the dark matter and stellar density (p) profiles, respectively, in the hydro- 
dynamic simulation. The thick black line corresponds to the dark matter density profile for the 
dark-matter-only simulation. The magenta line shows the velocity anisotropy, r\ (see caption to 
Figure 2 for the definition), profile for the dark matter (in the hydrodynamic simulation). 

center (0)- Two mechanisms contribute to the effect: first, the feedback flattens the central cusp, 
which reduces the efficiency of dynamical friction in the central regions ©; second, stellar 
feedback would have continued to heat the globular cluster orbits until stars stopped forming, 
around 200 Myr ago in Fornax ([9]). 

The distribution of velocities is isotropic within the core, and shows slight radial anisotropy 
outside the core (Figure 4), whereas the core remains isotropic throughout the evolution (Fig- 
ure 2). This behavior is inconsistent with a mechanism (1231) employing massive gas clouds, 
passively orbiting (not driven by feedback) near the galactic center, which flatten the dark mat- 
ter cusp via heating due to dynamical friction. It has been shown (|27l) that this would result 
in the development of significant tangential anisotropy within the core, which is not observed 
in our simulations. On the other hand, the gravitational resonance heating ((76|) naturally pro- 
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duces isotropic cores due to the fact that the feedback-driven bulk gas motions have random 
directions. 

These results also provide a natural explanation for the stellar population gradients seen in 
many early-type dwarfs (0|9]). In our simulations star formation is concentrated toward the 
galactic center. Over time, feedback gradually heats the population of stars, resulting in older 
(and more metal-poor) stars being kinematically warmer and having a larger spatial extent than 
younger (and more metal-rich) stellar populations. Hence we can reproduce, qualitatively, the 
age, metallicity, and velocity dispersion gradients observed in dwarf galaxies. 

Our simulations were stopped at z = 5 as continuing beyond this point would require a much 
larger computational box to correctly model the growth of larger structures and an infeasible 
increase in computation time. Furthermore, the impact of external ionizing radiation, ignored 
in our model, can become significant after z = 6.5. Nevertheless, we can reasonably infer the 
subsequent evolution of our model galaxy. If it is to become one of the early-type galaxies in 
the local universe (which are gas-poor), some mechanism will have to remove most or all of its 
interstellar medium. Some combination of a powerful star burst, increased metagalactic ionizing 
radiation, and ram-pressure stripping could result in the dwarf losing most of its gas (28). It is 
also likely that only a fraction of its star clusters will survive until the present time. As a result, 
our model galaxy would end up resembling a large dwarf spheroidal galaxy in the local universe: 
low stellar density; metal-poor with old stellar populations having pronounced radial population 
gradients; large stellar and dark matter cores (comparable in size and density to those in dwarf 
spheroidals); and perhaps a few globular clusters. In many respects, the galaxy would resemble 
the Fornax dwarf. 

Our non-cosmological modeling (16) suggested that stellar feedback can be directly respon- 
sible for the absence of dark matter cusps only in small galaxies, with total masses < 10 10 M : 
in larger galaxies the dark matter particle velocities become significantly larger than the veloc- 
ity of the random gas bulk motions, ~ 10 km s _1 . Our current, cosmological simulations are 
consistent with this result (the mass of our galaxy reaches 2 x 10 9 M by z = 5). Numerical 
simulations ([291) have suggested that a universal halo density profile (either cuspy or cored), 
once set is preserved through subsequent hierarchical evolution (which is consistent with the 
analytical result that the collisionless dark matter phase-space density can only decrease over 
time), implying that our mechanism may also lead to dark matter cores in larger galaxies. 

Our simulations indicate that the gravitational heating of matter resulting from feedback- 
powered bulk gas motions is a critical determinant of the properties of dwarf galaxies. Large 
dark matter cores are an unavoidable consequence of early star formation in dwarf galaxies. Our 
model indicates that, in primordial dwarf galaxies, globular clusters are formed in the most nat- 
ural place — near the center, where the gas pressure is highest — and are then pushed by feedback 
to much larger distances. This mechanism also ensures that globular clusters and unclustered 
stars have a comparable distribution, as is observed in early-type dwarfs (TJOl) . Additionally, the 
low stellar density and stellar population gradients observed in dwarf galaxies are also expected 
from the model. Finally, it is also worth noting that large cores have serious implications for 
direct searches of dark matter, as a flat core will produce a much weaker gamma ray annihilation 
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signal than a cusp. 
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Model 

We used GASOLINE (1571) , a parallel TreeSPH particle code that has individual particle timesteps 
In several respects our run parameters were similar to those used for the galaxy formation simu- 
lations run with GASOLINE reported in (S2). In what follows we will highlight the differences 
required to simulate small galaxies at high redshift. The code includes Compton cooling, atomic 
cooling based on collisional equilibrium and low-temperature cooling due to metals. The cos- 
mic ultraviolet background is assumed to be low at the redshifts we considered (S3) and was 
omitted. 

We model the low-temperature (< 10 4 K) radiative cooling of gas through fine structure 
and metastable lines of C, N, O, Fe, S, and Si following the prescription in (15?1) . To derive 
the cooling function, the authors assumed that the above elements are maintained in ionization 
equilibrium by locally produced cosmic rays, and that the cosmic ray ionization rate is scaled 
from the Galactic value by Z/Z Q , where Z/Z & is the metallicity of gas in solar units. We found 
that for temperatures T = 20 . . . 10 4 K their cooling function can be approximated very well by 
the following empirical expression: 

log(A/n|) = -24.81 + 2.928a; - 0.6982x 2 + log(Z/Z Q ), 

where x = log (log (log (T))). Here n u is the number density of hydrogen atoms (in cm -3 ), and 
A/n^ is in erg cm 3 s _1 units. 

The star formation algorithms we used were very similar to those that were extensively 
tested by Stinson et al. (1571) and Governato et al. (1521 . In order for gas to form stars it has to be 
Jeans unstable, to be colder than 15, 000 K, to have a minimum physical density of 10 atoms/cc 
and to have a minimum overdensity of 1,000. Star creation would then proceed at an average 
rate, in terms of solar masses per unit volume per unit time, of 0.05 times the gas density divided 
by the dynamical time. Star particles were created stochastically, as described in (15*51) , with a 
mass of 120 solar masses each. The probability that a star would be created in any million year 
interval was proportional to the average star formation rate of the parent gas particle multiplied 
by the time interval divided by the star mass. 

A Kroupa (S6) initial stellar mass function was used to determine the rate of supernovae per 
solar mass of stars. Each supernova was assumed to deposit -Esn = 0.4 x 10 51 ergs into the gas 
in net (allowing for some loss in the process). Our high mass resolution made it inappropriate 
to average the feedback effects due to a large population of stars when describing the feedback 
due to an individual star particle. Instead, we employed stochastic feedback by translating 
the time-averaged feedback into a probability that an individual supernova would occur in that 
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time. Thus the energy is injected in discrete supernova events. For larger star particle masses 
the energy output converges to the result with continuous feedback. Star particles select nearby 
gas particles in the manner described in (55) to receive feedback energy. However, the energy 
was added in a volume-weighted rather than mass-weighted manner which better respects the 
isotropic nature of the energy injection. 

Initial conditions 

In our simulations, we used the following values of the cosmological parameters (1571) : matter 
density VL m = 0.27, baryonic density £l b = 0.044, Hubble constant H = 71 km s _1 Mpc -1 , 
amplitude of fluctuations a 8 = 0.84, and spectral index n s = 0.93. We assumed that the 
universe is flat (fi m + tt^ = I, where £7a is the cosmological constant.) 

Constrained cosmological initial conditions were constructed using the package COSMICS 
(l5ffl) . We used two concentric spherical Gaussian ball constraints. The first constraint had 
a Gaussian scale length of R\ = 0.119 co-moving Mpc and an initial overdensity of 8\ = 
1.686D(0) / D(zcai\), where D(z) is the linear growth factor, and z co n = 5.92 is the targeted 
collapse redshift. This was designed to produce a halo with the mass of ~ 10 9 M by z ~ 6. 
After a few tests we realized that for our large box size (4 co-moving Mpc) one constraint was 
not enough to produce an isolated dwarf galaxy by z = 5. We chose a second constraint, with 
scale length of R 2 = 0.357 co-moving Mpc and initial overdensity of 82 = 5i[0.5(i?2/-Ri) 2 + 
O.5] 3 / 2 = 0.0894*1. (The value of 8 2 assumes an isolated density peak in an empty universe.) 

We used COSMICS to produce a cube populated with 1024 3 dark matter particles at zq = 
150. The central sphere, with radius of 0.4 co-moving Mpc and containing ~ 4.6 x 10 6 dark 
matter particles, was then populated with the same number of gas particles (each gas particle 
initially overlaying a dark particle and being assigned the same velocity). Following (1591) . the 
initial temperature of gas was set to 2.73(1 + z ) 2 /(l + 200) K= 310 K. (The above expres- 
sion assumes that the temperature of gas is the same as that of cosmic microwave background 
radiation until z = 200 and is set by adiabatic expansion thereafter.) 

To make the simulations computationally feasible, we lowered the resolution outside the 
central sphere by averaging coordinates and velocities for groups of adjacent particles. Between 
radii 0.4 and 0.8 co-moving Mpc, the number of particles was reduced by a factor of 2 3 , and 
in the remainder of the computational box by a factor of 8 3 . These low-resolution particles 
are needed to reproduce tidal torques from large scales on structures inside the central high- 
resolution sphere. Particles in the low resolution regions are collisionless and account for the 
mass of both dark matter and gas in those regions. 

The standard numerical code does not have the required physics (such as non-equilibrium 
molecular hydrogen chemistry and radiative transfer) to self-consistently describe star forma- 
tion in small progenitor galaxies with virial temperature T vir < 10 4 K. Proper modeling of these 
galaxies is important as they pre-enrich the intergalactic medium (out of which our model dwarf 
galaxy will be built) with trace amounts of heavy elements. We circumvented this limitation by 
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running the hydrodynamic simulations with a simplified prescription for star formation until 
z = 9.6 (before our model galaxy was assembled), and then switching to our default star for- 
mation recipe. The simplified portion of the simulations had no low-temperature (< 10 4 K) 
cooling due to heavy elements, and employed simple gas density (overdensity > 100 and den- 
sity > 1 cm -3 ) and temperature (T < 3 x 10 4 K) criteria to select star-forming particles. At 
z = 9.6, our model galaxy had a heavy element abundance of ~ 0.001 solar units, which is 
comparable to the lower stellar metallicity cutoff observed in nearby early-type dwarf galax- 
ies (15701) . (At the end of the simulations, the average metallicity of gas in the model galaxy 
reached 0.014 solar units.) Despite these significant simplifications, in the initial phase of the 
simulations we produce the correct abundance of heavy elements in the intergalactic medium. 
We also explored a range of heavy elements abundances at z = 9.6, by scaling the metallicity 
of all gas particles by a constant factor in the range £ = 0.2 to 5 (see the next section). Af- 
ter z = 9.6, we ran the simulations with the full physics (low-temperature cooling and Jeans 
instability criterion). 

Numerical tests 

To estimate the smallest resolved radius at the center of our model galaxy, we carried out addi- 
tional dark-matter-only simulations which had 8 times lower mass resolution than our fiducial 
model. The initial conditions for this run were generated from the initial conditions for the high 
resolution run by binning together groups of 8 adjacent particles (reducing the effective resolu- 
tion by a factor of two in each linear dimension). We chose the gravitational softening length 
for the new run to be twice larger (e = 24 pc) than for the fiducial run. 

Figure SI shows the final radial density profiles for both high and low resolutions runs. 
One can see that both profiles converge around radius ~ 80 pc (or 3.3e), implying that our 
high resolution model should resolve radii down to ~ 40 pc. The actual resolution is probably 
somewhat smaller in the fiducial run. At the radius where the low and high resolution models 
start diverging, the low resolution run exhibits both a break in the density profile (becomes 
shallower) and the one-sigma spread becomes significantly wider (see Figure SI). For the high 
resolution run, the same behavior is observed at a radius of 25 pc, or 2e. Thus we are confident 
that in the fiducial models we resolve central radii to at least 40 pc, and perhaps to 25 pc. This is 
sufficient for our purposes, as stellar feedback significantly perturbs the distribution of galactic 
gas on scales > 100 pc. 

Next we explored model sensitivity to the three important model parameters: E$n (the en- 
ergy deposited into the interstellar medium by a single supernova); the low-temperature cooling 
factor, k (see the description below); and the initial metallicity factor, £. (The fiducial values 
are E SN = 0.4 x 10 51 ergs, k — 1, and £ = 1.) It would be computationally prohibitive to 
run full-scale simulations for different values of the above parameters. Instead, we ran a se- 
quence of small-box simulations. At z = 9.6, we cut out a sphere centered on our model galaxy 
with radius 6.1 kpc (three virial radii at that epoch). We then ran eight different cosmological 
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Figure SI: Radial density profiles for the model galaxy for dark-matter-only simulations at 
z = 5. Solid blue and red lines show the dark matter density profiles averaged over the last 50 
snapshots for the fiducial model and the model with 8 times lower mass resolution, respectively. 
The dotted blue and red lines show the corresponding one-sigma deviations from the averaged 
profiles. The vertical solid green line is at 80 pc, the radius where the low and high resolution 
models converge. The vertical dashed green line corresponds to the estimated smallest resolved 
radius (40 pc) for the high resolution run. 

simulations, with the extracted sphere placed in an empty universe. Two of the runs had identi- 
cal parameters to the fiducial full-box simulations (a dark-matter-only run and a hydrodynamic 
run). Six other small-box simulations were run with two different values for each of the three 
free parameters, with the rest of the parameters kept at their fiducial values. 

In Figures S2-4 we show the evolution of the central dark matter phase space density, F, 
measured in the same way as in the full-box simulations (see the main text of the paper). Specif- 
ically, it was measured for the central 5000 dark matter particles, typically located within the 
central ~ 200 pc. This quantity is relatively insensitive to the variable adiabatic compression 
of dark matter caused by baryons (gas and stars). Black and green lines correspond to fidu- 
cial values of the parameters for the dark-matter-only and hydrodynamic small-box simulations 
respectively (these lines are repeated in all three figures). 
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Figure S2: Evolution of the central dark matter phase- space density, F, for small-box simula- 
tions with different values of E SN . Black and green lines correspond to the dark- matter-only 
and hydrodynamic fiducial models (E$n = 0.4 x 10 51 ergs). Blue and red lines correspond to 
the models with four times lower (0.1 x 10 51 ergs) and four times larger (1.6 x 10 51 ergs) values 
of the parameter E SN , respectively. 
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Figure S3: Evolution of the central dark matter phase space density, F, for small-box simu- 
lations with different values of k. Black and green lines correspond to the dark-matter-only 
and hydrodynamic fiducial models (k = 1). Blue and red lines correspond to the models with 
k = 0.2 and 5, respectively. 
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Figure S4: Evolution of the central dark matter phase space density, F, for small-box simu- 
lations with different values of £. Black and green lines correspond to the dark-matter-only 
and hydrodynamic fiducial models (£ = 1). Blue and red lines correspond to the models with 
£ = 0.2 and 5, respectively. 
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Figure S5: Evolution of the mass-to-light ratio within the central 200 pc, (M/ L) , for small-box 
simulations with different values of E^n- Blue, green, and red lines correspond to the models 
with Esn = (0.1, 0.4, 1.6) x 10 51 ergs, respectively. 

A comparison of Figures S2-4 with the bottom panel of Figure 2 from the main text, shows: 
(a) for the dark-matter-only models, in the small-box simulations the quantity F behaves sim- 
ilarly to the case of the full-box simulations (converges to an approximately constant value); 
and (b) gravitational dark matter heating due to stellar feedback manifests itself in both large- 
and small-box simulations, but the magnitude of the effect is smaller in the smaller box case. 
The latter result is not unexpected, as in the small-box run the amount of dark matter and gas 
available for the build-up of the model galaxy is much smaller than in the full-box simulations, 
resulting in much less energetic star formation after the initial strong star burst around z ~ 7.7. 

Reducing the energy input from supernovae by a factor of four does not affect appreciably 
gravitational heating of dark matter (see Figure S2). Four times larger energy input results in 
somewhat stronger effect. The effect is clearly seen for the whole range of Esn tested. 

Figure S5 demonstrates that we produce quite dark galaxies, with the central value of the 
total mass-to-light ratio, (M/L) , reaching 2, 5, and 9 by z = 5 for the models with E SN = 
(0.1, 0.4, 1.6) x 10 51 ergs, respectively. (In calculating (M/L) we adopted the mass-to-light 
ratio of 1 for stars (I5ii|) , and ignored gas, as current dwarf spheroidal galaxies have very little 
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gas.) These values are close to the observed ones; e.g., Leo I and Fornax have (M/L) = 3 and 
5, respectively (1572I) . 

As discussed in § 1, the low-temperature (T < 10 4 K) radiative cooling function we use in 
our model is based on several assumptions. We investigated the sensitivity of the assumptions 
by running small-box simulations with the original cooling function multiplied by a constant 
factor, k, equal to 0.2, 1 , and 5. The impact of stellar feedback on the central phase space density 
of dark matter is quite insensitive to the efficiency of low-temperature cooling as shown in 
Figure S4. We emphasize, however, that to properly resolve star formation and stellar feedback 
in dwarf galaxies, some form of low-temperature cooling has to be present in the model. In our 
runs with no radiative cooling for T < 10 4 K, the shortest Jeans length for gas particles inside 
the galaxy was ~ 10 kpc - much larger that the size of the galaxy. As a result, the galactic 
disk was absolutely gravitationally stable, and any star formation prescription produced stars 
at random locations inside the disk. In this case stars are not formed in clusters and stellar 
feedback does not produce large-scale gas bulk motions (instead heating the interstellar gas 
almost uniformly). Not surprisingly, in such models the dark matter experiences essentially no 
gravitational heating due to the feedback. 

The last parameter we explored was the initial abundance of heavy elements inside the 
model galaxy. The fiducial gas metallicity at z = 9.6 was ~ 0.001 solar units (see § 1). We 
ran two small-box models with the metallicity of each gas particle multiplied by a factor £ 
equal to 0.2 and 5. Figure S5 shows that both the lower and higher metallicity cases resulted 
in a somewhat stronger gravitational heating effect. This may be due to either a significant 
non-linear dependence of the strength of the effect on gas metallicity, or (more likely) on the 
stochastic nature of our star formation and stellar feedback prescriptions. Again, for the range 
of £ tested, gravitational heating due to stellar feedback is significant. 

To summarize, the small-box simulations demonstrated that our main result — that stellar 
feedback results in significant gravitational heating of dark matter at the centers of primordial 
dwarf galaxies — is largely insensitive to the values of the model parameters £sn> and £ when 
they vary within reasonable ranges. We believe that this is because small galaxies are effectively 
"pressure cookers" in which star formation and feedback are confined to the central parts of the 
halos. A combination of model parameters that results in more energy being pumped into the 
interstellar gas will result in individual star bursts that have a stronger effect on dark matter, but 
the interval between star bursts becomes longer (it takes longer for hotter gas to cool down), so 
the overall cumulative effect is qualitatively similar. 

It follows from our numerical experiments that to accurately model formation and evolution 
of a dwarf galaxy in a cosmological context, the following conditions have to be met: (a) The 
model has to include a prescription for low-temperature (< 10 4 K) radiative gas cooling due 
to heavy elements. This makes the simulations significantly (by at least a factor of 10) more 
expensive, but is absolutely necessary if one wants to to resolve the critical star formation and 
feedback processes, (b) In SPH codes, the gas mass resolution has to be high enough (a few 
hundred M , or better) to properly resolve the extremely low density gas inside supernovae 
remnants: it is the thermal pressure of this hot (> 10 6 K) gas that drives the large-scale bulk 
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motions of the interstellar gas. 
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